#ifndef DRAINAGE_HPP
#define DRAINAGE_HPP

#include "meshTypes_fwd.hpp"
#include "solverBaseTypes_fwd.hpp"
#include "sweSolver.hpp"
#include <experimental/memory>

namespace hms
{

template<typename MeshType>
struct Drainage {
private:
	std::experimental::observer_ptr<const SWESolver<MeshType>> m_solver {nullptr};
	std::vector<Vector2s> m_manholeCoords;

public:
	Drainage() = default;

	Drainage( const SWESolver<MeshType>& solverArg ){
		solver( solverArg );
	}

	const SWESolver<MeshType>& solver() const {
		assert( m_solver && "Dereferencing nullptr!" );
		return *m_solver;
	}
	void solver( const SWESolver<MeshType>& solverArg ){
		m_solver = std::experimental::make_observer( &solverArg );
	}

	void operator() (ArrayXXs& vals){
		/* currently, the implementation relies on uniform cell sizes to find
		 * the cells in which manholes are located. Therefore, it won't work
		 * for other types */
		assert( (std::is_same_v<MeshType, UniMesh>) );
		if constexpr ( std::is_same_v<MeshType, UniMesh> ){
			const UniMesh& mesh { solver().mesh().derived() };
			for ( const Vector2s& coords : m_manholeCoords ){
				Array2i cellIndices {
					( ( coords - mesh.origin() ).array() / mesh.diff() )
					.template cast<Index>()
				};
				[[maybe_unused]]
				UniCell cell { mesh.cell(cellIndices[0], cellIndices[1]) };
				vals = 0;
			}
		}
	}
};

template<typename MeshType>
void setDrainageCallable( SWESolver<MeshType>& solver ){
	using SFName = typename SWESolver<MeshType>::SourceFieldName;
	if ( solver.useSourceField( SFName::Drainage ) ){
		auto drainageFunctor { Drainage(solver) };
		SourceField<NonUniform, MeshType> drainageField {
			solver.mesh(), std::move(drainageFunctor), 1
		};
		drainageField.storage() = 0;
		solver.sourceField(SFName::Drainage) = std::move(drainageField);
	}
}


} // namespace hms

#endif